Working with HDF5 files in R

Link to tutorial page

Set up and view data file

library(rhdf5)
library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
##  Path to GDAL shared files: C:/Users/sjgraves/Documents/R/win-library/3.3/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: C:/Users/sjgraves/Documents/R/win-library/3.3/rgdal/proj
##  Linking to sp version: 1.2-3

Get information from H5 file

# set file path
# use ../ to move up a directory to the data folder
h5file <- "../NEONdata/D17-California/TEAK/2013/spectrometer/reflectance/Subset3NIS1_20130614_100459_atmcor.h5"

# view h5 structure
h5ls(h5file)
##    group                                 name       otype  dclass
## 0      /                     ATCOR_Input_File H5I_DATASET  STRING
## 1      /                 ATCOR_Processing_Log H5I_DATASET  STRING
## 2      /                Aerosol Optical Depth H5I_DATASET INTEGER
## 3      /                               Aspect H5I_DATASET   FLOAT
## 4      /                          Cast Shadow H5I_DATASET INTEGER
## 5      / Dark Dense Vegetation Classification H5I_DATASET INTEGER
## 6      /                 Haze-Cloud-Water Map H5I_DATASET INTEGER
## 7      /                  Illumination Factor H5I_DATASET INTEGER
## 8      /                          Path Length H5I_DATASET   FLOAT
## 9      /                   Processing Version H5I_DATASET  STRING
## 10     /                          Reflectance H5I_DATASET INTEGER
## 11     /                Shadow_Processing_Log H5I_DATASET  STRING
## 12     /                      Sky View Factor H5I_DATASET INTEGER
## 13     /               Skyview_Processing_Log H5I_DATASET  STRING
## 14     /                                Slope H5I_DATASET   FLOAT
## 15     /          Slope_Aspect_Processing_Log H5I_DATASET  STRING
## 16     /                  Solar Azimuth Angle H5I_DATASET   FLOAT
## 17     /                   Solar Zenith Angle H5I_DATASET   FLOAT
## 18     /                    Surface Elevation H5I_DATASET   FLOAT
## 19     /                 Visibility Index Map H5I_DATASET INTEGER
## 20     /                   Water Vapor Column H5I_DATASET INTEGER
## 21     /             coordinate system string H5I_DATASET  STRING
## 22     /                       flightAltitude H5I_DATASET   FLOAT
## 23     /                        flightHeading H5I_DATASET   FLOAT
## 24     /                           flightTime H5I_DATASET   FLOAT
## 25     /                                 fwhm H5I_DATASET   FLOAT
## 26     /                             map info H5I_DATASET  STRING
## 27     /              to-sensor azimuth angle H5I_DATASET   FLOAT
## 28     /               to-sensor zenith angle H5I_DATASET   FLOAT
## 29     /                           wavelength H5I_DATASET   FLOAT
##                dim
## 0                1
## 1                1
## 2    544 x 578 x 1
## 3    544 x 578 x 1
## 4    544 x 578 x 1
## 5    544 x 578 x 1
## 6    544 x 578 x 1
## 7    544 x 578 x 1
## 8    544 x 578 x 1
## 9                1
## 10 544 x 578 x 426
## 11               1
## 12   544 x 578 x 1
## 13               1
## 14   544 x 578 x 1
## 15               1
## 16           1 x 1
## 17           1 x 1
## 18   544 x 578 x 1
## 19   544 x 578 x 1
## 20   544 x 578 x 1
## 21               1
## 22         5732053
## 23         5732053
## 24         5732053
## 25         426 x 1
## 26               1
## 27   544 x 578 x 1
## 28   544 x 578 x 1
## 29         426 x 1
# import spatial information
mapInfo <- h5read(h5file,
                  "map info",
                  read.attributes = T)
# set reflectance info
reflInfo <- h5readAttributes(h5file,
                             "Reflectance")

# set scale factor and no data value objects
scaleFactor <- reflInfo$`Scale Factor`

noDataValue <- as.numeric(reflInfo$`data ignore value`)

Use functions to automatically extract this information

NEON wrote functions to do many data extractions manually. Information is found on the website

rm(list = ls())

source("neon_hdf5_functions.R")

f <- "../NEONdata/D17-California/TEAK/2013/spectrometer/reflectance/Subset3NIS1_20130614_100459_atmcor.h5"

Apply functions to TEAK data

# define the CRS definition by EPSG code
epsg <- 32611

# set wavelengths
wavelengths <- h5read(f,"wavelength")
### final Code ####
# H5close()

# find the dimensions of the data to help determine the slice range
# returns cols, rows, wavelengths
dims <- get_data_dims(fileName = f)

# open band, return cleaned and scaled raster
band <- open_band(fileName=f,
                  bandNum = 56,
                  epsg=epsg)

# plot data
plot(band,
     main="Raster for Lower Teakettle - B56")

# extract 3 bands
# create  alist of the bands
bands <- list(58, 34, 19)

# use lapply to run the band function across all three of the bands
rgb_rast <- lapply(bands, open_band,
                   fileName=f,
                   epsg=epsg)

# create a raster stack from the output
rgb_rast <- stack(rgb_rast)

# plot the output, use a linear stretch to make it look nice
plotRGB(rgb_rast,
        stretch='lin')

# CIR create  alist of the bands
bands <- c(90, 34, 19)

CIRStack <- create_stack(f, 
                         bands, 
                         epsg)
plot_stack(CIRStack,
           title="Color Infrared (CIR) Image")

# create a list of the bands
bands <- list(152,90,58)
aStack <- create_stack(f, bands, epsg)
plot_stack(aStack,
           title="another combo")

Save data

# export as a GeoTIFF
writeRaster(CIRStack,
            file="../outputs/TEAK/TEAK_CIR_image.tif",
            format="GTiff",
            overwrite=TRUE)